Copyright © 1973 American Telephone and Telegraph Company 

The Bell System Technical Journal 

Vol. 52, No. 7, September, 1973 

Printed in U.S.A. 



Distribution of Ha n /n, a n Randomly 
Equal to ±1 
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(Manuscript received February 6, 1973) 

When a% } a 2 , ■ • • are independent random variables, each equal to -±.1 
with ■probability \, the sum £" a„/n is a random variable whose distribu- 
tion is difficult to determine theoretically. This sum is of interest in the 
study of intersymbol interference in digital communication systems. Here 
the distribution of the sum is computed by numerical integration and the 
results tabulated. Asymptotic expressions are given for the tails of the 
distribution. 

I. INTRODUCTION 

The distribution of the random variable 

00 

X = £ ajn, (1) 

n— 1 

where a h a 2 , • • • are independent random variables equal to +1 or 
— 1 with probability ^, is of some interest in the study of intersymbol 
interference in a digital communication system. For example, the sum 
of two independent expressions of the form (1) occurs when the pulse 
train £-«o a n sin (t — nir)/(t — nir), a n randomly equal to ±1, is 
sampled at regularly spaced instants which are slightly out of step with 
the zeros of sin t. The theory of random variables of type (1) (in par- 
ticular with a„/3" in place of a n /n) has been studied by a number of 
investigators. A survey of the field has been made recently by Hill and 
Blanco. 1 Here we evaluate the distribution of x numerically and give 
expressions for its behavior when x is large. Questions of continuity 
and convergence are put aside. 

Since the distribution is even about x = 0, only values for x ^ 
need be considered. From the characteristic function 

°0 

f(u) = avg [exp (ixu)~\ = IJ [cos (u/n)~] (2) 

n-l 
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we get an expression for the probability density p(x) of x: 

1 roo 

p(x) = - / cos {xu)f(u)du, (3) 

7T Jo 

Prob (a; > *0 = o / 7. — /( w ) dw - ( 4 ) 

Z 7T ./0 W 

The values of p(rc) and Prob (x > zi) shown in Table I were obtained 
by evaluating these integrals by the trapezoidal rule 23 which works well 
for (3) and (4). 

The asymptotic expressions (8) and (18) for p(x) follow from a 
saddle point analysis of (3). Both p(x) and Prob (x > Xi) decrease 
rapidly when x (or Xi) > 3, the decrease being dominated by the factor 
exp [— exp(x — A)2 where A = 1.39 

The rapid decrease of p(x) is interesting because the divergence of 
£ 1/n might lead one to expect that p(x) would decrease slowly as 
x— >°°. Instead, p(x) actually decreases much faster than a Gaussian 
probability density. I am indebted to a referee for the observation 
that the second, fourth, and sixth moments of p(x) are, respectively, 
?r 2 /6, ll7r 4 /180, and 233tt 6 /7560. 

The reader may wonder why as many as six decimal places are given 
in Table I. There are several reasons. One is that the cost was low. 
About 3 seconds were required by a Honeywell 6000 Processor to 
compute the values shown in Table I, and about 40 terms were required 
in each trapezoidal sum. This illustrates the fact (apparently not well 
known) that when integrals like (3) and (4) are to be evaluated 
numerically, the trapezoidal rule often performs better than most of 
the other conventional quadrature methods (better than Simpson's 
rule, for instance). 2 The six-figure accuracy is also used to gain an idea 
of the values of x for which the asymptotic expansion (18) for p(x) 
begins to be valid. This degree of accuracy also shows that p(0) is 
equal to 0.249 994 ■■ ■ and not to i, as might be inferred from a four- 
figure tabulation. 

II. TRAPEZOIDAL RULE CALCULATION 

Preliminary computations showed that | f(u) \ < 10~ 10 when u > 15. 
Furthermore, it was found that (3) and (4) could be evaluated to 
within the desired accuracy by using a trapezoidal-rule spacing of 
Aw = h = 0.4. In line with these values the trapezoidal sum was 
truncated at the 40th term (15/0.4 ^ 40). 
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Table I — Values of p(x) and Prob (x > Xi) 



x or xi 


P(x) 


Prob (x > Xi) 


x or X\ 


p(x) 


Prob (x > xi) 


0.0 


0.249 994 


0.500 000 


2.0 


0.125 000 


0.056 599 


0.2 


0.249 970 


0.450 003 


2.2 


0.091 768 


0.034 949 


0.4 


0.249 802 


0.400 021 


2.4 


0.061 647 


0.019 683 


0.6 


0.249 073 


0.350 118 


2.6 


0.037 148 


0.009 912 


0.8 


0.246 778 


0.300 494 


2.8 


0.019 592 


0.004 357 


1.0 


0.241 222 


0.251 623 


3.0 


0.008 777 


0.001 623 


1.2 


0.230 408 


0.204 357 


3.2 


0.003 222 


0.000 494 


1.4 


0.212 852 


0.159 912 


3.4 


0.000 927 


0.000 118 


1.6 


0.188 353 


0.119 683 


3.6 


0.000 198 


0.000 021 


1.8 


0.158 232 


0.084 949 


3.8 


0.000 029 


0.000 003 


2.0 


0.125 000 


0.056 599 


4.0 


0.000 003 


0.000 000 



The infinite product (2) for f(u) was computed by using 

f(u) = ( II [cos (u/n)l) exp £ In cos (u/n) , (5) 

\ n-l / .n-JV J 

where N is a large number such that u/N <3C 1 for all values of u used 
in the computation (0 ^ u ^ 16). The product Ili' -1 m (5) was 
computed by straightforward multiplication. The sum in (5) was 
computed by setting g{n) = In [cos (u/n)2 in the Euler-Maclaurin 
sum formula : 

t 9(n) = f g(t)dt + \g(N) ~ fj gW(N) 

n=N JN & 4 1 



-f> 3) (A0 



B, 



(2k) I 



g<**-»(N) +R k . (6) 



Here the B's denote Bernoulli's numbers, B 2 = 1/6, 2? 4 = — 1/30, 
B 6 = 1/42 [we stopped at B 6 in our use of (6) J, g (k) (N) denotes the 
value of (d/dt) k g(t) at t = N, and the remainder R k is the integral of 
(2 * +1) (<) times the Bernoulli "polynomial" of degree 2k + 1 and 
period 1 (see pages 520-540 of Ref. 4). 

The integral in (6) can be evaluated to within the desired accuracy 
by setting u/t = y, dt = — udy/y 2 , expanding In (cos y) in powers of 
y with the help of 



/v rv r v 3 

tan vdv = \v -f — + 



2v 5 17t; 7 
15 + 315 ~ t " 



dv, 
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and integrating termwise: 

roo fu/N 

/ g(t)dt = u \ y~ 2 In (cos y)dy 

J N JO 

~ 2iV T 36# 3 T 225JV 5 T 17640W 7 T ' w 

Expressions for the higher derivatives of g(N) in (6) can be obtained by 
differentiating the series in (7) with respect to N. 

In using (6) we stopped at k = 3 and neglected R 3 . 

Three separate trapezoidal-rule evaluations of p{x) and Prob (x > x{) 
were made using eqs. (3) to (7) with h = 0.4, N = 201; h = 0.38, 
N - 201; and h = 0.36, N = 301, respectively. Here h is the spacing 
used in the trapezoidal-rule evaluations of (3) and (4). The three sets 
of computed values differed only in the 7th or 8th decimal places, i.e., 
all agreed with the values shown in the table. The values 201 and 301 
of N are so large that terms beyond k = 3 in (6) and those shown in 
(7) are not needed. To check the computations, the integrals of p(x) 
amd x 2 p(x) from x = to x = « (the upper limit used was actually 
x = 5) were computed by the trapezoidal rule with a spacing of 
Ax = 0.1. The trapezoidal values agreed with the known values, 
respectively $ and tt 2 /12, to within 6 significant figures or better. 

III. DISCUSSION OF TABLE I 

Table I shows that p(x) remains nearly equal to p(0) = 0.249994 for 
£ £ £ 1, passes through p(2) = 0.125000 (is it exactly £?), and then 
decreases rapidly. The question as to whether p(2) is exactly | remains 
unanswered, but p(0) = 0.249994 does not seem to be an erroneous 
calculation of \. For if p 2 (u) is the probability density of u = L2 ajn, 
then 

p(x) = \p%ix - 1) + hp*(x + 1). 

Setting x equal to and 2 and combining the results give a result I 
owe to J. E. Mazo, 

p(2) = Jp(0) + ip»(3) > Jp(0). 
Furthermore, replacing £p 2 (3) by p(4) — §?>2(5) gives 
p(2) = M0) + p(4) - *p.(5) 
0.125000 = 0.124997 + 0.000003 - M(5) 

which is satisfied by the tabulated values when p 2 (5) is assumed to 
be negligibly small. 
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IV. ASYMPTOTIC EXPRESSIONS FOR LARGE X 

We shall show that the rapid decrease of p(x) for x > 3 is described 
by 

2/o = cxp O - 2 7 + In (tt/4)] = exp [re - 1.39599 • • •], ^ 

where y denotes Euler's constant, 0.577215 • • • . Integrating (8) gives 

Prob (x > xx) - erfc (yj) ~ p(xi)/y , (9) 

where y is computed from the second of equations (8) with X\ in place 
of .r. Bounds for the distribution involving exponential functions of 
e x have been obtained by L. A. Shepp in unpublished work. 
To obtain (8) we rewrite the integral (3) for p(x) as 

1 f K 
P(«) = 9- / e iiu f(u)du. (10) 

As is often the case for such integrals, the asymptotic value of p(x) is 
given by the contribution of a saddle point, u = iy , lying far out on 
the positive imaginary w-axis (the path of integration being deformed 
so as to pass through the saddle point). For u = iy the integrand in 
(10) becomes 

exp[-.r?/+ <p(y)], tp(y) = £ In [cosh (y/n)]. (11) 

When y is large we can show that 

v(v) = V In y - y + Ay + i In 2 + r(y), 
A = 2 7 - In (tt/4) = 1.39599 • • •, (12) 

where r(y) has roughly the same magnitude as exp (~2y). 

We first outline the derivation of the expression (12) for <p(y), and 
then apply (12) to obtain the asymptotic expression (8) for p(x). 

The derivation of (12) is based upon the Euler-Maclaurin sum 
formula (6) with X = 1 and g(t) = In [cosh {y/t)~\. The integral in the 
sum formula is 

J In [cosh (y/t)~\dt 

= y / v~ 2 In (cosh v)dv 
Jo 

= - In (cosh y) + t/(ln y) tanh y - y f (In v) sech 2 vdv 

Jo 

= ln2~y-\-y\ny-yf (In v) sech 2 vdv + 0[ye- 2 " In 2/], (13) 
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where we have integrated by parts twice. The last integral in (13) has 
the value 

f* (In v) sech 2 vdv = - y + In (x/4) (14) 

Jo 
which can be obtained by (i) replacing In v in (14) by v"- 1 , (ii) differ- 
entiating the known value (formula 3.527-3, page 352 of Ref. 5) of the 
resulting integral with respect to n, and (Hi) setting p - 1. The 
derivative of g(t) with respect to t, g {1) (t), in the sum formula (6) is 

g ii)(t) = - yt-i tanh (y/t) 

= _ yt -i(i - 2e- 2 «" + 2e~ i « lt ), 

where the exponential terms become negligible when y becomes large 
and t = 1. In general, for I = 0, 1, 2, • • • , 0< 2I+1) (1) is equal to 
— (21 + 1) ! y plus negligible terms. Therefore, the right side of the sum 
formula (6) is the integral plus 

J(y-Jn2)+§ 1 y + ^V+---+l* (15) 

plus terms which are negligible when y is large. The sum of the co- 
efficients of y in (15) is known to be equal to y (page 529 of Ref. 4). 
Hence (15) is equal to yy - \ In 2. Addition of (13) and (15) and use 
of (14) gives the expression (12) for <p(y). 

Next we use the expression for tp(y), with the small term r(y) 
neglected, to obtain the asymptotic form of p(x). The saddle point of 
interest occurs at u = iyo where y is the zero of the derivative of the 
exponent in (11). The exponent is -xy + <p(y), and y is the zero of 

-x + <p'(y) =- x+\ny + A. 

Thus 2/0 = exp (x — A). This y Q is the same as the y appearing in the 
asymptotic expression for p(x) stated in (8). The exponent itself has 
the value -xy + <p(yo) = - 2/0 + 5 In 2 at y . By making use of 
p'Hy) = \/y and the higher derivatives of tp(y), the exponent can be 
expanded in a Taylor series about y . From this expansion it follows 
that near u the integrand in the integral (10) for p(x) can be written as 

«p [-»+ 3 ha + •£#%]• (16) 

where z = u - u . Setting (16) in (10), changing the variable of 
integration from u to t = z/yo = (u — u )/yo, and assuming that p(x) 
is given asymptotically (as x and y tend to <») by the contribution of 
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the saddle point at u tl give 

p(x) - (2T)- I 2*y e-» f exp [ y.j£ fj^r^ ] dr. (17) 

Here the nominal path of integration is the real r-axis. The classical 
saddle point asymptotic expansion obtained from (17) is 



p(x) ~ (Wir) 1 *-" 



1 + 2^-nHs + ---! "" 



The coefficients of the powers of l/y in the series can be determined by 
a general procedure described in Appendix D, page 1999, of Ref. 6. 
The asymptotic expression stated in (8) is the leading term in (18). 
An idea of the accuracy of the asymptotic expressions can be obtained 
by considering the case x = 3. For x = 3, yo is 4.973 and Table I gives 
the "exact" value p(3) = 0.008777. The asymptotic value of p(3) 
obtained from (8) is 0.008710, the first two terms in (18) give 0.008783, 
and the first three terms give 0.008776. Table I gives the "exact" 
value Prob (x > 3) = 0.001623 and eq. (9), namely Prob (x > 3) 
-erfc (j/i) = erfc (2.230), gives 0.001612. 
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